/*
 * =====================================================================================
 *
 *       Filename:  ncmatrix.h
 *
 *    Description:
 *
 *        Version:  1.0
 *        Created:  2009年06月02日 15时10分48秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Ying Wang (WY), ywang@nlpr.ia.ac.cn
 *        Company:  Institute of Automation, Chinese Academy of Sciences
 *
 * =====================================================================================
 */

#ifndef NCMATRIX_H
#define NCMATRIX_H

#include <iostream>
#include <sstream>
#include "mathutils.h"
#include "ncvector.h"
template<typename T> class NCmatrix;

/**
 *
 *
 */
template<typename T>
class NCmatrix
{
	template<typename T1>
		friend std::ostream& operator << (std::ostream& os, const NCmatrix<T1>& mat);
	template<typename T1>
		friend std::istream& operator >> (std::ostream& is, NCmatrix<T1>& mat);

public:
	NCmatrix( );
	NCmatrix( int row, int column );
	NCmatrix( int row, int column, T a );
	NCmatrix( int row, int column, const T *a );
	NCmatrix( const NCmatrix<T>& mat );
	NCmatrix<T>& operator = ( const NCmatrix<T>& mat );
public:
	~NCmatrix();

public:
	inline NCmatrix<T> transpose () const;
public:
	inline NCmatrix<T>& operator += ( const T &a );
	inline NCmatrix<T>&  operator *= (const T &a);

	inline NCmatrix<T>& operator += ( const NCmatrix<T> &mat );
	inline NCmatrix<T>& operator -= ( const NCmatrix<T> &mat );
	inline NCmatrix<T>  operator +  ( const NCmatrix<T> &mat ) const;
	inline NCmatrix<T>  operator -  ( const NCmatrix<T> &mat ) const;

	inline NCmatrix<T>  operator *  (const T &a) const;
	inline NCmatrix<T>  operator *  ( const NCmatrix<T> &mat ) const;
	inline NCvector<T>  operator *  ( const NCvector<T> &vec ) const;

public:
	inline T* operator [] (int index_row);
	inline const T* operator [] (int index_row) const;
	inline const T operator () ( int index_row, int index_column ) const;
	inline T& operator () ( int index_row, int index_column ) ;
	inline T** getData(){return data;}
	inline const T** getData() const {return data;}
	inline void setValue( int irow, int jcolumn, const T &value);
public:
	inline int row() const
	{
		return m_row;
	}
	inline int column() const
	{
		return m_column;
	}
    //template<typename T>
	inline NCvector<T> rowvector(int row) const;
	inline NCvector<T> columnvector(int column) const;
	inline void resize(int row, int column);

	inline NCmatrix<T> submat(int irow,int jcolumn, int row, int column) const;
	inline NCmatrix<T> slicemat(const NCvector<bool> &slice) const;
public:
	inline NCmatrix<T>* doubleSizeMatrix() const;
	inline NCmatrix<T>* doubleSizeMatrix2() const;
private:
	void freeData();
private:
	int m_row;
	int m_column;
	T **data;

public:
	inline static NCmatrix<T> eye( int n );
};

template<class T>
NCmatrix<T> :: NCmatrix()
	: m_row( 0 )
	, m_column( 0 )
	, data( NULL )
{

}

template<typename T>
NCmatrix<T> :: NCmatrix( int row, int column )
	: m_row( row )
	, m_column( column )
	, data( m_row>0 ? new T*[m_row] : NULL)
{
	int numele = m_row * m_column;
	if(data != NULL)
	{
		data[0] = numele>0 ? new T[numele] : NULL;
		for(int i=1; i<m_row; i++)
		{
			data[i] = data[i-1] + m_column;
		}
	}

}

template<typename T>
NCmatrix<T> :: NCmatrix( int row, int column, T a)
	: m_row( row )
	, m_column( column )
	, data( m_row>0 ? new T*[m_row] : NULL)
{
	int  i,numele = m_row * m_column;
	if( data != NULL)
	{
		data[0] = numele>0 ? new T[numele] : NULL;
		for( i=0; i<numele; i++)
		{
			data[0][i] = a;
		}
		for( i=1; i<m_row; i++ )
		{
			data[i] = data[i-1] + m_column;
		}
	}
}

template<typename T>
NCmatrix<T> :: NCmatrix( int row, int column, const T *a)
	: m_row( row )
	, m_column( column )
	, data( m_row>0 ? new T*[m_row] : NULL)
{
	int i, numele = m_row * m_column;
	if ( data != NULL )
	{
		data[0] = numele>0 ? new T[numele] : NULL;
		for ( i=0; i<numele; i++)
		{
			data[0][i] = a[i];
		}
		for( i=1; i<m_row; i++)
		{
			data[i] = data[i-1] + m_column;
		}
	}
}
template<typename T>
NCmatrix<T> :: NCmatrix( const NCmatrix<T>& mat)
	: m_row( mat.m_row )
	, m_column( mat.m_column )
	, data( m_row>0 ? new T*[m_row] : NULL)
{
	int i,numele = m_row * m_column;
	if ( data != NULL )
	{
		data[0] = numele>0 ? new T[numele] : NULL;
		for ( i=0; i<numele; i++ )
		{
			data[0][i] = mat.data[0][i];
		}
		for ( i=1; i<m_row; i++ )
		{
			data[i] = data[i-1] + m_column;
		}
	}
}

template<typename T>
NCmatrix<T>& NCmatrix<T> :: operator = ( const NCmatrix<T>& mat )
{
	if( this != &mat )
	{
		int i, numele;
		if( m_row == mat.m_row && m_column == mat.m_column && data != NULL)
		{
			numele = m_row * m_column;
			for(i=0; i<numele; i++)
			{
				data[0][i] = mat.data[0][i];
			}
		}
		else
		{
			freeData();
			m_row = mat.m_row;
			m_column = mat.m_column;
			numele = m_row *m_column;
			data = m_row>0 ? new T*[m_row] : NULL;
			if( data != NULL)
			{
				data[0] = numele>0 ? new T[numele] : NULL;
				for( i=0; i<numele; i++ )
				{
					data[0][i] = mat.data[0][i];
				}
				for( i=1; i<m_row; i++ )
				{
					data[i] = data[i-1] + m_column;
				}
			}

		}
	}
	return *this;
}

template<typename T>
NCmatrix<T> :: ~NCmatrix( )
{
	freeData();
}
template<typename T>
inline NCmatrix<T>& NCmatrix<T> :: operator += ( const T &a )
{
	int i, numele = m_row * m_column;
	for ( i=0; i<numele; i++ )
	{
		data[0][i] += a;
	}
	return *this;
}

template<typename T>
inline NCmatrix<T> NCmatrix<T> :: transpose () const
{
	NCmatrix<T> result(m_column,m_row);
	for ( int i=0; i<m_row; i++ )
	{
		for ( int j=0; j<m_column;j++ )
		{
			result.data[j][i] = data[i][j];
		}
	}
	return result;
}

template<typename T>
inline NCmatrix<T>& NCmatrix<T> :: operator *= ( const T &a )
{
	int i, numele = m_row * m_column;
	for ( i=0; i<numele; i++)
	{
		data[0][i] *= a;
	}
	return *this;
}

template<typename T>
inline NCmatrix<T>& NCmatrix<T> :: operator += ( const NCmatrix<T> &mat )
{
	if ( m_row == mat.m_row && m_column == mat.m_column )
	{
		int i, numele = m_row * m_column;
		for ( i=0; i<numele; i++ )
		{
			data[0][i] += mat.data[0][i];
		}
		return *this;
	}
	else
	{
		std::cerr << "The size of the two matrices must equal!" << std::endl;
		throw( "The size of the two matrices must equal!" );
	}
}

template<typename T>
inline NCmatrix<T>& NCmatrix<T> :: operator -= ( const NCmatrix<T> &mat )
{
	if ( m_row == mat.m_row && m_column == mat.m_column )
	{
		int i, numele = m_row * m_column;
		for ( i=0; i<numele; i++)
		{
			data[0][i] -= mat.data[0][i];
		}
		return *this;
	}
	else
	{
		std::cerr << "The size of the two matrices must equal!" << std::endl;
		throw( "The size of the two matrices must equal!" );
	}
}

template<typename T>
inline NCmatrix<T> NCmatrix<T> :: operator + ( const NCmatrix<T> &mat ) const
{
	if ( m_row == mat.m_row && m_column == mat.m_column )
	{
		NCmatrix result(m_row,m_column);
		int i, numele = m_row * m_column;
		for ( i=0; i<numele; i++ )
		{
			result.data[0][i] = data[0][i] + mat.data[0][i];
		}
		return result;
	}
	else
	{
		std::cerr << "The size of the two matrices must equal!" << std::endl;
		throw( "The size of the two matrices must equal!" );

	}
}

template<typename T>
inline NCmatrix<T> NCmatrix<T> :: operator - ( const NCmatrix<T> &mat ) const
{
	if ( m_row == mat.m_row && m_column == mat.m_column )
	{
		NCmatrix result(m_row,m_column);
		int i, numele = m_row * m_column;
		for ( i=0; i<numele; i++ )
		{
			result.data[0][i] = data[0][i] - mat.data[0][i];
		}
		return result;
	}
	else
	{
		std::cerr << "The size of the two matrices must equal!" << std::endl;
		throw( "The size of the two matrices must equal!" );

	}
}

template<typename T>
inline NCmatrix<T> NCmatrix<T> :: operator * (const T &a) const
{
	int i, numele = m_row * m_column;
	NCmatrix<T> result(m_row,m_column);
	for( i=0; i<numele; i++ )
	{
		result.data[0][i] = data[0][i]*a;
	}
	return result;
}

template<typename T>
inline NCmatrix<T> NCmatrix<T> :: operator * ( const NCmatrix<T> &mat) const
{
	if ( m_column == mat.m_row )
	{
		NCmatrix result(m_row,mat.m_column);
		int i,j,k;
		double tmp;
		for ( i=0; i<m_row; i++ )
		{
			for ( j=0; j<result.m_column; j++ )
			{
				tmp = 0;
				for ( k=0; k<m_column; k++ )
				{
					tmp += data[i][k] * mat.data[k][j];
				}

				result.data[i][j] = tmp;
			}
		}
		return result;
	}
	else
	{
		std::cerr << "These two matrices cannot multiple!" << std::endl;
		throw( "These two matrices cannot multiple!" );
	}
}

template<typename T>
inline NCvector<T> NCmatrix<T> :: operator * ( const NCvector<T> &vec) const
{
	if( m_column == vec.size() )
	{
		NCvector<T> result(m_row);
		int i,j;
		double tmp;
		for( i=0; i<m_row; i++ )
		{
			for( tmp=0.,j=0; j<m_column; j++ )
				tmp += data[i][j] * vec[j];
			result[i] = tmp;
		}
		
		return result;
	}
	else
	{
		std::cerr << "Matrice cannot multiple the vector" << std::endl;
		throw( "Matrice cannot multiple the vector");

	}
}

template<typename T>
inline T* NCmatrix<T> :: operator [] ( int index_row )
{
	if ( index_row >= 0 && index_row < m_row )
	{
		return data[index_row];
	}
	else
	{
		std::cerr << "Index out of range!" << std::endl;
		throw( "Index out of range!" );
	}
}

template<typename T>
inline const T* NCmatrix<T> :: operator [] (int index_row) const
{
	if ( index_row >= 0 && index_row < m_row )
	{
		return data[index_row];
	}
	else
	{
		std::cerr << "Index out of range!" << std::endl;
		throw( "Index out of range!" );
	}
}
template<typename T>
inline const T NCmatrix<T> :: operator () ( int index_row , int index_column ) const
{
	if( index_row >=0 && index_row < m_row )
	{
		return data[index_row][index_column];
	}
	else
	{
		std::cerr << "operator(), index out of range!" <<std::endl;
	}
}
template<typename T>
inline T& NCmatrix<T> :: operator () ( int index_row, int index_column ) 
{
	if( index_row >=0 && index_row < m_row )
	{
		return data[index_row][index_column];
	}
	else
	{
		std::cerr << "operator(), index out of range!" <<std::endl;
	}

}

template<typename T>
inline void NCmatrix<T> :: setValue ( int irow, int jcolumn, const T &value)
{
	assert( irow >=0 && irow < m_row );
	assert( jcolumn >=0 && jcolumn < m_column );
	data[irow][jcolumn] = value;
}

template<typename T>
void NCmatrix<T> :: freeData( )
{
	if ( data != NULL )
	{
	    if(data[0] != NULL)
            delete [](data[0]);
		delete []data;
		m_row = 0;
		m_column = 0;
	}
}

template<typename T>
inline NCvector<T> NCmatrix<T> :: columnvector(int column) const
{
	if(column>=0 && column<m_column)
	{
		NCvector<T> result(m_row);
		for(int i=0; i<m_row; i++ )
			result[i] = data[i][column];
		return result;
	}
	else
	{
		std::cerr<< "matrix::columnvector(int column), index out of range"<<std::endl;
		throw("matrix::columnvector(int column), index out of range");
	}
}

template<typename T>
inline NCvector<T> NCmatrix<T> :: rowvector(int row) const
{
    if(row>=0 && row<m_row)
    {
        NCvector<T> result(m_column);
        for(int j=0;j<m_column;j++)
            result[j] = data[row][j];
        return result;
    }
    else
    {
        std::cerr<<"matrix::rowvector(int row), index out of range"<<std::endl;
        throw("matrix::rowvector(int row), index out of range");
    }
}

template<typename T>
inline NCmatrix<T> NCmatrix<T> :: submat(int irow, int jcolumn, int row, int column) const
{
	if( irow+row > m_row || jcolumn+column > m_column )
	{
		std::cerr<<"submat, out of range"<<std::endl;
		throw("submat,out of range");
	}
	NCmatrix<T> result(row,column);
	for(int i=0; i<row; i++ )
	{
		for(int j=0; j<column; j++ )
		{
			result[i][j] = data[i+irow][j+jcolumn];
		}
	}
	return result;
}

template<typename T>
inline NCmatrix<T> NCmatrix<T> :: slicemat(const NCvector<bool> &slice) const
{
	if(slice.size() != m_row)
		std::cerr<<"NCmatrix<T> :: slicemat, the size of slice are not consistent with matrix"<<std::endl;
	int i,j,numtrue = 0,islice=0;
	for( i=0; i<m_row; i++ )
	{
		if(slice[i])
			numtrue++;
	}
	NCmatrix<T> result(numtrue,m_column);
	for( i=0; i<m_row; i++ )
	{
		if(slice[i])
		{
			for( j=0; j<m_column; j++ )
			{
				result[islice][j] = data[i][j];
			}
			islice++;
		}
	}
	return result;
}
template<typename T>
inline NCmatrix<T>* NCmatrix<T> :: doubleSizeMatrix() const
{
	int row = 2 * m_row, column = 2 * m_column;
	NCmatrix<T>* mat = new NCmatrix<T>( row, column );
	for( int i=0; i<row; i++ )
	{
		for( int j=0; j<column; j++ )
			mat->setValue( i, j, data[i/2][j/2]);

	}
	return mat;
}

template<typename T>
inline NCmatrix<T>* NCmatrix<T> :: doubleSizeMatrix2() const
{
	int row = 2 * m_row, column = 2 * m_column;
	NCmatrix<T>* mat = new NCmatrix<T>( row, column );
	int irow, jcolumn;
	for( irow=0; irow<row; irow+=2 )
	{
		for ( jcolumn=0; jcolumn<column; jcolumn+=2 )
		{
			mat->setValue(irow,jcolumn,data[irow/2][jcolumn/2]);
		}
	}
	for( irow=0;irow<row; irow++ )
		mat->setValue(irow,column-1,data[irow/2][(column-1)/2]);
	for( jcolumn=0; jcolumn<column; jcolumn++ )
		mat->setValue(row-1,jcolumn,data[(row-1)/2][jcolumn/2]);

	/*
	 *  A B C
	 *  E F G
	 *  H I J
	 */
	for( irow=1; irow<row-1; irow += 2 )
	{
		for ( jcolumn=0; jcolumn<column; jcolumn += 2 )
		{
			mat->setValue(irow, jcolumn, (data[irow/2][jcolumn/2]+ data[irow/2+1][jcolumn/2])/2.0);
		}
	}

	for( irow=0; irow<row; irow+=2 )
	{
		for( jcolumn=1; jcolumn<column-1; jcolumn+=2 )
		{
			mat->setValue(irow,jcolumn,(data[irow/2][jcolumn/2]+data[irow/2][jcolumn/2+1])/2.0);
		}
	}
	for( irow=1; irow<row-1; irow+=2 )
	{
		for( jcolumn=1; jcolumn<column-1; jcolumn+=2 )
		{
			mat->setValue(irow,jcolumn,(data[irow/2][jcolumn/2]+data[irow/2][jcolumn/2+1]
						   +data[irow/2+1][jcolumn/2]+data[irow/2+1][jcolumn/2+1])/4.0);
		}
	}
	return mat;	
}

template<typename T>
inline void NCmatrix<T> :: resize(int row, int column)
{
	if(m_row != row || m_column != column)
	{
		int i, numele=row*column;
		freeData();
		m_row = row;
		m_column = column;
		data = m_row>0 ? new T*[m_row] : NULL;
		if( data != NULL)
		{
			data[0] = numele>0 ? new T[numele] : NULL;

			for( i=1; i<row; i++ )
				data[i] = data[i-1] + m_column;
		}


	}
}

template<typename T1>
std::ostream& operator << (std::ostream &os, const NCmatrix<T1> &mat)
{
	for(int i=0; i<mat.m_row; i++)
	{
		for(int j=0; j<mat.m_column;j++)
		{
			os << mat.data[i][j] <<" ";
		}

		os<<std::endl;
	}
	return os;
}

template<typename T1>
std::istream& operator >> (std::istream &is, NCmatrix<T1> &mat)
{
	std::string strline;
	getline(is,strline);
	std::stringstream stream(strline);
	int i,j,row=0, column =0;
	T1 tmp;
	while(stream>>tmp)
		column ++ ;
	stream.clear();
	row++;
	while(getline(is,strline))
		row++;
	mat.resize(row,column);
	is.clear();
	is.seekg(0,std::ios::beg);

	for( i=0; i<row; i++ )
	{
		for( j=0; j<column; j++ )
		{
			is>>tmp;
			mat[i][j] = tmp;
		}
	}

	return is;

}

template<typename T>
inline NCmatrix<T> NCmatrix<T>::eye( int n )
{
	T a=0;
	NCmatrix<T> result(n,n,a);
	for ( int i=0; i<n; i++)
	{
		result.data[i][i] = 1;
	}
	return result;
}
#endif
